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A  numerical  modeling  approach  was  developed  to  simulate  the  propagation  of 
shear  and  longitudinal  sound  waves  in  arbitrary,  dense  dispersions  of  spherical 
particles.  The  scattering  interactions  were  modeled  with  vector  multipole 
functions  and  boundary  condition  solutions  for  each  particle.  Multiple  scattering 
was  simulated  by  translating  the  scattered  wave  fields  from  one  particle  to  another 
with  the  use  of  translational  addition  theorems,  summing  the  multiple-scattering 
contributions,  and  recalculating  the  scattering  using  an  iterative  method.  The 
theory  and  initial  results  for  the  model  are  presented,  including  an  integral 
derivation  for  the  translational  addition  theorems.  The  model  can  simulate  3D 
material  microstructures  with  a  variety  of  particle  size  distributions,  compositions, 
and  volume  fractions.  To  test  the  model,  spectra  and  wave  field  patterns  were 
generated  from  both  ordered  and  disordered  microstructures  containing  up  to 
several  hundred  particles.  The  model  predicts  wave  propagation  phenomena  such 
as  refractive  focusing,  mode  conversion,  and  band  gap  phenomena.  The 
convergence  of  the  iterations  ranges  from  excellent  to  fair,  and  is  dependent  on 
the  field  (longitudinal  or  shear),  particle  configuration,  and  elastic  wave 
frequency.  The  model  is  currently  limited  by  the  computation  of  sufficiently  high 
multipole  order  for  the  simulation  of  dense  particle  dispersions. 

PACS  numbers:  43.20.Gp,  43.35.Cg 

I.  INTRODUCTION 

The  propagation  of  elastic  waves  through  particles  randomly  dispersed  in  a  liquid  or  solid  is 
relevant  to  many  applications,  including  the  nondestructive  evaluation  of  particulate  composites, 
the  design  of  advanced  acoustic  materials,  the  process  monitoring  of  suspensions,  the  ultrasonic 
examination  of  living  tissues,  and  seismic  models  for  various  geologic  media  (sediments,  rocks, 
etc.).  Modeling  elastic  waves  in  such  systems  is  difficult,  however,  due  to  multiple  scattering 
between  particles,  mode  conversion  from  scattering  at  particle  surfaces,  lack  of  periodicity  in  the 
particle  configuration,  and  heterogeneity  of  particle  sizes  and  compositions.  Non-empirical, 
physics-based  models  are  therefore  desired  to  simulate  elastic  wave  propagation  in  dense, 
random  dispersions  of  particles  and  predict  the  resultant  macroscopic  wave  propagation 
properties. 

This  paper  presents  the  theory  and  results  of  a  numerical  model  developed  to  address 
these  goals.  The  model  uses  vector  multipole  functions,  boundary  condition  solutions, 
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translational  addition  theorems,  and  direct  iteration  to  solve  the  multiple-scattering 
problem  for  an  arbitrary  distribution  of  spheres  in  a  solid  matrix.  The  model  can 
calculate  either  the  spatial  patterns  or  frequency  spectra  of  elastic  wave  fields  for  a 
simulated  random  particle  configuration.  Macroscopic  properties  such  as  attenuation  and 
velocity  can  therefore  be  predicted  for  a  wide  variety  of  particulate  systems. 

Multiple  scattering  in  particulate  systems  is  not  solvable  in  an  exact,  closed  analytical  form 
for  an  arbitrary  number  and  arrangement  of  particles.  Scattering  models  must  therefore  use 
either  statistical  formulations  or  numerical  simulation.  The  simulation  of  elastic  waves  in 
complex  particle  dispersions  is  computationally  exhausting,  however,  for  many  modeling 
approaches.  Finite  methods  are  the  most  widely  used  numerical  approaches  for  modeling  fields 
in  physical  systems,  but  are  often  limited  to  2D  representations  or  the  use  of  periodic  cells  for  the 
simulation  of  particle  dispersions.1"  '  These  simplifications  reduce  the  number  of  calculations, 
but  may  also  exclude  effects  arising  from  a  fully  random  and  3D  particle  configuration. 

Multipole  approaches  can  be  more  efficient  than  finite  methods  for  modeling  systems  with 
spherical,  spheroidal,  or  cylindrical  particles  since  a  multipole  field  contains  significantly  more 
information  than  a  field  value  at  a  grid  point.  Several  multipole  techniques  have  been  developed 
for  multiple  scattering,  including  transfer  matrix  (T  matrix)  methods,  fast  multipole  methods,  and 
direct  iterative  methods  (as  opposed  to  iterative  solution  methods  for  matrix  approaches). 

The  direct  iterative  method  has  been  applied  to  electromagnetic  scattering  in  small 
collections  of  spherical  particles,  where  it  is  also  known  as  the  order-of-scattering  method.4"6 
Direct  iterative  methods  have  yet  to  be  applied,  however,  to  fully  random  3D  particle  systems 
with  large  numbers  of  particles  (>100  particles),  high  particle  concentrations  (>20%),  and  mixed 
particle  sizes  and  properties.  Additionally,  multipole-based  acoustic  models  for  particulate 
systems  have  been  limited  to  either  longitudinal  wave  propagation  or  periodic  lattices  of 
particles.  "  Little  work  has  been  reported  on  the  modeling  of  full  elastic  wave  scattering  in 
random  particulate  systems,  with  longitudinal  and  shear  waves  in  both  particles  and  matrix. 

There  are  a  number  of  reasons  to  justify  the  use  of  a  direct  iterative  approach  over  a  T  matrix, 
fast  multipole,  or  unit-cell  approach.  The  T  matrix  approach  solves  the  multiple-scattering 
problem  as  a  linear  system  of  equations.  In  contrast,  the  direct  iterative  approach  solves  for 
multiple  scattering  by  iteratively  recomputing  the  multiple-scattering  contributions  until 
convergence  of  the  elastic  wave  fields.  Direct  iterative  methods  are  therefore  computationally 
simpler  and  avoid  the  pitfalls  of  numerically  inverting  large  matrices  for  large  numbers  of 
particles.  Often,  such  matrices  are  solved  using  iterative  methods,  which  could  be  argued  as 
placing  the  direct  iterative  approach  on  the  same  level  as  the  T  matrix  approach.  However,  the 
direct  iterative  approach  sidesteps  the  need  to  explicitly  formulate  and  solve  the  multiple 
scattering  as  a  problem  in  linear  algebra. 

An  additional  advantage  of  the  iterative  approach  is  its  correspondence  to  the  physically 
intuitive  order-of-scattering  concept,  where  each  iteration  step  represents  a  successive  order  of 
scattering.  For  example,  the  initial  fields  before  iteration  correspond  to  the  zero-order 
contributions  to  multiple  scattering  (single-scattering  approximation),  the  first  iteration 
corresponds  to  the  first-order  contributions,  the  second  iteration  corresponds  to  the  second-order 
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contributions,  and  so  on.  The  ability  to  separately  compute  and  inspect  each  successive 
scattering  order  provides  a  deeper  understanding  of  the  scattering  process  and  of  the  effects  of 
multiple  scattering.  It  also  provides  a  useful  diagnostic  of  the  model’s  performance  and 
accuracy,  and  allows  certain  interactions  to  be  excluded  for  computational  efficiency  (e.g., 
interactions  between  distant  particles). 

In  comparison  to  fast  multipole  methods,  the  iterative  approach  treats  the  particle 
configurations  in  a  more  straightforward  manner  without  resorting  to  segregation  of  the  particles 
into  groups  with  short-range  (particle-to-particle)  and  long-range  (group-to-group)  interactions. 
The  fast  multipole  method  reduces  the  computations  associated  with  such  simulations  by 
imposing  a  hierarchical  scattering  process  that  assumes  material  homogeneity  at  a  sufficiently 
large  scale  above  the  particle  level.10  This  assumption  limits  the  applicability  of  the  approach  to 
macroscopically  homogeneous  composites  with  little  spatial  variation  in  the  microstructure.  The 
iterative  multipole  method,  however,  imposes  no  such  restrictions  on  the  material  microstructure, 
and  is  additionally  numerically  simpler. 

For  ordered  particle  systems,  unit  cell  methods  employ  periodic  boundary  conditions  to 
model  an  infinite  lattice,  and  this  method  has  also  been  applied  to  random  ensembles.  The 
limitations  of  this  method  include  the  difficulty  of  modeling  finite  lattices  (<103  particles), 
defects,  and  various  levels  of  disorder.  The  unit  cell  approach  may  therefore  be  insufficient  for 
modeling  acoustic  band  gap  phenomena  for  point  and  line  defects,  lattice  heterostructures,  and 
small  acoustic  band  gap  (phononic)  crystals.  When  applied  to  random  configurations  (i.e.,  the 
unit  cell  is  comprised  of  a  random  particle  arrangement),  the  method  introduces  long-range  order 
that  may  bias  the  model’s  results. 

This  paper  presents  the  theory  and  initial  results  for  an  iterative  multipole  model  developed 
for  elastic  wave  propagation  in  particulate  systems.  The  particles  are  modeled  as  spheres 
embedded  in  a  solid  matrix,  permitting  the  use  of  spherical  wave  functions  for  the  elastic  wave 
fields.  The  model  computes  the  propagation  of  waves  by  using  single  particle  scattering 
solutions,  addition  theorems  to  translate  the  scattered  fields  from  one  particle  to  another,  and 
iteration  to  compute  the  scattering  resulting  from  the  multiply  scattered  fields.  Simulations  were 
performed  on  a  personal  computer  for  particle  systems  of  up  to  several  hundred  particles  and 
particle  volume  fractions  up  to  50%.  Initial  results  are  presented  for  example  simulations,  and 
the  capabilities  and  deficiencies  of  the  current  model  are  summarized.  Methods  for  improving 
both  the  accuracy  and  efficiency  of  the  model  are  discussed. 

II.  THEORY 

A.  Vector  multipole  functions 

Spherical  wave  expansions  form  the  theoretical  basis  for  the  iterative  multipole  method,  and 
describe  elastic  waves  derived  from  the  Navier  equation  for  linear,  homogeneous  materials: 

p—~Y  =  (■%  +  2//)V(V  ■  u)  -  //V  x  (V  x  u)  .  (1) 

dt 
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Eq.  (1)  can  be  recast  as  two  separate  vector  Helmholtz  equations  by  splitting  the  displacement 
vector  u  into  longitudinal  and  transverse  parts  and  imposing  an  hannonic  time  dependency  of  the 
form  e~,cot : 


(V2  +k2L)wL  =0, 

(2) 

(V2+k2)u,=0. 

(3) 

Solutions  to  the  vector  Helmholtz  equation  in  spherical  coordinates  are  the  vector  spherical 
wave  functions.  The  vector  spherical  wave  functions  in  this  work  were  constructed  from 
spherical  radial  functions  and  vector  spherical  harmonics  denoted  as  zn  ( kr )  and  Y'mn  ( 6 ,  (p) 
respectively.  These  functions  are  also  referred  to  as  vector  multipole  functions.  The  spherical 
radial  functions  are  comprised  of  the  spherical  Bessel  functions  j  n  ( kr )  for  standing  waves 

inside  the  particles,  spherical  Hankel  functions  of  the  first  kind  h(nl>  (kr)  for  outward-propagating 
scattered  waves,  and  spherical  Hankel  functions  of  the  second  kind  h(n2)  (kr)  for  inward- 
propagating  incident  waves. 

The  vector  spherical  harmonics  used  in  this  work  are  eigenfunctions  of  the  quantum  angular 
momentum  operators  J2 ,  Jz ,  L2 ,  and  S2 ,  and  are  also  known  as  pure-orbital  vector  spherical 
harmonics.11'13  Using  helicity  basis  vectors  and  scalar  spherical  harmonics  they  can  be  defined 
as 


'm+ Ie-1 


Y L(0,<P)  =  CX+u,-lYlJ 

+  C"mXoYl,meO  + 


(4) 


The  C-symbols  are  the  Clebsch-Gordan  coefficients  commonly  used  in  quantum  mechanics  and 
derived  from  integrals  involving  three  spherical  harmonics.  The  notation  used  is  from 
Varshalovich  et  al.  Vector  functions  are  necessary  for  representing  shear  elastic  fields,  and  the 
vector  spherical  harmonics  provide  a  compact  notation  for  otherwise  awkward  combinations  of 
scalar  spherical  harmonics. 

Since  they  are  solutions  to  the  Helmholtz  equation,  the  vector  multipole  functions  can  be 
generated  from  a  scalar  potential  ®  and  vector  potential  *F  defined  as  following: 


® = Z  Z  z<  ikr]Y.~  (0.  <P)  (5) 


f=IZ".wY:w?)  («) 

72=0  m=—n 


The  vector  multipole  functions  are  derived  by  applying  the  gradient  and  curl  operations  to  the 
potentials  in  the  following  manner: 
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U  =  -V<D 
k 


(V) 


v  =  IvxT 

k 


(8) 


W  =  — rV  x(VxT) 
k 


(9) 


The  resulting  vector  multipole  functions  are  complete  and  orthogonal  multipole  expansions  for  n 
=  0  to  oo  and  m  =  -n  to  +n,  and  are  the  following: 


n—0  m=-n 


v 


2/2  +  1 


Y1  _i_  1 

+tez-(Ar)Y-1(^^ 


(10) 


v  =  XI 

«=0  m——n 


n  +  1 
2/?  +  l 


-2, 


2/7  +  1 


(id 


w  =  XZ^(^)Yl(^^) 


«=0  m——n 


(12) 


U,  V,  and  W  are  the  longitudinal,  electric,  and  magnetic  vector  multipole  functions, 
respectively.  They  are  similar  to  other  widely  used  wave  function  and  multipole  formulations, 
specifically  the  L,  N,  and  M  vector  spherical  wave  functions  of  Stratton14;  the  A(r;L),  A (/•;£), 
and  A(r;M)  vector  multipoles  of  Greiner  and  Maruhn,  Rose,  and  Edmonds1145,16;  and  the  X 
vector  spherical  harmonic  of  Jackson.17  The  U,  V,  and  W  functions,  however,  are  specifically 
consistent  with  the  definitions  in  Eqs.  (5-9)  as  opposed  to  other  vector  multipole  definitions,  and 
provide  a  more  compact  notation  than  Stratton’s  vector  spherical  wave  functions. 

U  corresponds  to  the  longitudinal  displacement  field.  V  and  W  are  transverse  fields  that  can 
correspond  to  either  elastic  shear  fields  or  electromagnetic  fields.  For  elastic  waves,  V 
corresponds  to  the  shear-electric  displacement  field  (shear  field  with  electric-type  vector 
multipole  function),  and  W  corresponds  to  the  shear-magnetic  displacement  field  (shear  field 
with  magnetic-type  vector  multipole  function).  The  horizontal  and  vertical  conventions  normally 
used  to  differentiate  the  two  types  of  shear  fields  are  not  used  here  due  to  the  ambiguity  of 
distinguishing  vertical  and  horizontal  orientations  in  a  spherical  coordinate  system. 
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The  initial  waves  incident  on  the  particles  are  modeled  as  plane  waves.  The  spherical 
expansion  for  the  longitudinal  plane  wave  is  derived  from  the  gradient  of  the  partial  wave 
expansion  for  scalar  fields: 


ikz 


1 

ik 


1V 

ik 


oo 

Y.i'‘j4rt2^7T)j,,(kr)Y,JIW,<P) 


=  Z  >'"'V4*(2'>  +  1)U„.0 

n= 0 


(13) 


At  least  three  different  partial  wave  expansions  have  been  presented  in  the  literature  for 
transverse  vector  planes  waves,  including  those  by  Stratton,14  Jackson,17  and  Greiner  and 
Maruhn. 1 1  Extensive  mathematical  analysis  and  numerical  testing,  however,  has  shown  that  the 
expansions  of  Greiner  and  Maruhn  present  the  correct  formulation  and  consistently  converge. 
These  expansions  were  reformulated  for  the  U,  V,  and  W  functions,  and  for  shear  waves 
propagating  in  the  z  direction  and  polarized  in  the  x  and  y  directions.  These  partial  wave 
expansions  are  the  following: 


e  e 


ikz 


eA 


ikz 


ZVA(2»  +  l)[w„.„  +V,M  +  W„,_,  -V,„],  (14) 

n 

(-;)£;■  A(  2»  +  l)[w„,  +  V„,  -  W„_,  +V,_,].  (15) 

n 


B.  Boundary  condition  solutions 

Iterative  simulation  of  multiple  scattering  in  a  dispersion  of  spherical  particles  first  requires 
scattering  solutions  for  the  individual  spheres.  For  scattering  from  a  single  solid  sphere  in  a  solid 
matrix,  there  will  be  an  incoming  incident  field,  a  refracted  internal  field,  and  an  outgoing 
scattered  field  for  each  of  the  U,  V,  and  W  wave  fields.  Each  of  these  wave  field  components 
will  also  have  an  associated  amplitude  coefficient.  Fig.  1  shows  the  relationship  between  each  of 
the  field  components  and  the  coefficients. 

Single-sphere  scattering  has  been  solved  numerous  times  in  the  literature,  but  often  not  in  the 
most  general  form.  These  solutions  mostly  model  the  incident  fields  as  plane  waves  with  only  a 
longitudinal  or  shear  component.  "  Although  these  conditions  are  often  sufficient  for  dilute 
suspensions,  the  incident  fields  for  each  particle  in  the  iterative  multipole  approach  will  not  in 
general  be  plane  waves,  but  a  combination  of  the  initial  plane  waves  and  the  scattered  waves 
from  other  particles.  Also,  the  incident  fields  will  be  a  combination  of  both  longitudinal  and 
shear  waves.  Single-sphere  scattering  solutions  were  therefore  required  for  arbitrary  incident 
fields  comprised  of  both  longitudinal  and  shear  waves. 


Approved  for  public  release;  distribution  unlimited. 


Iterative  simulation  of  scattering  6 


Longitudinal 
(white  arrows) 


nm  /v. 


E 


nm 


Shear-electric 
(grey  arrows) 


Shear-magnetic 
(black  arrows) 


G 


nm 


nm 


i 


nm 


Incident 

Waves 


Refracted 

(interior) 

Waves 


Scattered 

Waves 


FIG.  1 .  Diagram  of  incident,  refracted,  and  scattered  elastic  waves  for  single-particle  scattering, 
with  associated  amplitude  coefficients  for  the  (n,m)  multipole  moment. 

Given  an  arbitrary  incident  field,  the  amplitude  coefficients  of  the  refracted  and  scattered 
fields  are  found  by  solving  the  boundary  conditions  on  the  surface  of  the  sphere.  The  boundary 
conditions  provide  a  set  of  six  linear  equations  for  the  six  unknown  coefficients.  Three  of  these 
equations  are  obtained  from  continuity  of  the  displacements,  and  the  other  three  are  derived  from 
continuity  of  the  stresses,  where  u  is  the  vector  displacement  and  a  is  the  stress  tensor: 


^incident  ^ scattered  _ ^ refracted 


(16) 


_ incident  ,  _ scattered  _ refracted 

<j  +  <J  =  cr 


(17) 


Orthogonality  conditions  for  the  vector  spherical  harmonics  are  used  to  eliminate  the  angular 
dependence  from  the  boundary  conditions.  The  resulting  equations  have  the  following  matrix 
form  which  relate  the  six  unknown  coefficients  to  the  three  known  coefficients: 


V\  (j) 

-7iW 

72  0) 

-7  20) 

0 

0  " 

( D  ) 

^ NM  1 

f  dNM  7l(g')  +  ^NM  72  is ) 

7«0) 

“7«W 

7  7  O') 

-7  70) 

0 

0 

^ NM 

d^NM  7f,  is)  BnM  dl  is) 

74  O’) 

“74W 

75  O') 

~  7,5  0) 

0 

0 

77 

^  NM 

d^MddS^^^NMdiiS) 

790) 

-79W 

7,o  0) 

-7,oO) 

0 

0 

H nm 

d  Kim  V9(s)  +  B  nm  7 1 0  O  ) 

0 

0 

0 

0 

73  0) 

-  73  0) 

F 

1  NM 

^  NM^I  3  (s  ) 

l  0 

0 

0 

0 

7s  0) 

-7gO)y 

V  I  NM  J 

v  Ow7s(g)  y 

The  77-symbols  are  functions  of  the  spherical  radial  functions,  multipole  order  N,  longitudinal 
wave  vector  Ay,  shear  wave  vector  ks,  and  the  sphere  radius  a.  The  j,  h,  and  g  in  the  /7-functions 
refer  to  the  type  of  radial  function  used  in  ?;  [  j  =  j  n  ( kr ) ,  h  =  hff  ( kr ) ,  and  g  =  hf]  (kr)].  The  77- 
functions  are  given  in  the  appendix,  with  the  appropriate  radial  function  substituting  for  z.  Since 
the  solution  matrix  separates  into  a  4x4  matrix  and  2x2  matrix,  the  shear-magnetic  fields  are 
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decoupled  from  the  longitudinal  and  shear-electric  fields  and  do  not  participate  in  mode 
conversion.  The  scattering  amplitudes  for  each  sphere  are  obtained  in  the  scattering  model  by 
solving  the  2x2  and  4x4  matrices. 

C.  Translation  of  scattered  fields 

After  computing  the  wave  fields  due  to  single-sphere  scattering  for  all  of  the  particles  in  the 
simulation,  the  scattered  waves  from  each  particle  are  then  propagated  to  the  other  particles  to 
become  part  of  their  incident  wave  fields.  New  scattered  waves  are  then  recomputed  with  use  of 
the  above  boundary  conditions,  and  the  process  is  repeated.  This  procedure  iteratively  updates 
the  incident  fields  with  the  multiple-scattering  contributions  until  the  amplitudes  of  the  scattered 
waves  converge.  Since  the  scattered  wave  fields  for  each  particle  are  in  the  original  particle’s 
coordinate  system,  the  vector  multipole  functions  for  each  particle  must  be  transfonned  into  the 
coordinate  systems  of  the  other  particles  in  order  to  update  the  incident  fields.  This  is 
accomplished  with  the  use  of  special  spherical  wave  expansions  known  as  translational  addition 
theorems.  The  transfonnation  of  the  vector  multipole  functions  from  one  coordinate  system  to 
another  are  expressed  as  follows: 


co  v 


u  =y  y  onm  u' 

^  nm  zLsVfj.  ^  v/d  ’ 

(19) 
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v=0  /u=-v 

(21) 

Fig.  2  shows  the  geometric  relationship  for  the  field  translations  between  two  spherical 
particles.  The  untranslated  fields  UBm,  \nm  ,  and  W„m  in  Eqs.  (19)-(21)  represent  waves 

scattered  from  the  original  or  “transmitting”  sphere  a,  and  are  in  sphere  a’s  coordinate  system. 
The  translated  fields  U'  ,  V  ,  and  represent  the  scattered  waves  incident  on  a  second  or 

“receiving”  sphere  [!„  and  are  transformed  to  sphere  P’s  coordinate  system.  The  global  position 
vectors  for  the  two  spheres  are  R(/  and  Rp.  The  position  of  sphere  a  with  respect  to  sphere  p  is 
therefore  R«p=  R„  -  Rp.  The  translation  vector  ra  is  in  sphere  a’s  local  coordinate  system.  The 
translational  addition  theorems  therefore  express  the  vector  multipole  functions  in  the  a 
coordinate  system  (ra,  6a,  and  <pa)  as  expansions  of  multipole  functions  in  the  P  coordinate 
system  (rp,  dp,  and  cpp).  These  translation  coefficients  are  only  valid  on  the  surface  of  sphere  P 
(rp=  ap,  where  ap  is  the  radius  of  sphere  P).  Note  that  the  center  of  sphere  a  must  lie  outside  of 
sphere  P  (ap  <  Rap ).  This  restricts  the  use  of  the  addition  theorems  to  spheres  external  to  each 
other  and  non-overlapping.  Although  addition  theorems  can  also  be  derived  for  spheres 
embedded  within  larger  spheres  (ap  >  Rap),  this  work  will  only  concern  itself  with  the  more 
commonly  used  external  forms. 
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Sphere  (3 


FIG.  2.  Position  (R„  and  Rp)  and  translation  (ra)  vectors  for  addition  theorems  with  respect  to 
local  and  global  coordinates  for  spheres  a  and  p. 

The  scalar  addition  theorem  is  sufficient  for  translating  the  longitudinal  field  since  the 
longitudinal  vector  multipole  function  arises  from  a  scalar  potential  [Eq.  (7)].  The  scalar 
addition  theorem  has  been  published  extensively,24'27  and  provides  the  Q"™  translation 
coefficients: 
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n,m 

v,ju 
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(22) 


The  vector  addition  theorems  are  required  to  translate  the  shear  fields  since  the  transverse 
vector  multipole  functions  arise  from  a  vector  potential.  Cruzan  was  one  of  the  earliest  to  derive 
expressions  for  the  Snv”  and  T™  translation  coefficients  by  directly  translating  the  coordinates  in 

the  vector  spherical  wave  functions  and  applying  various  identities  and  relationships  to  arrive  at 

28 

an  analytical  solution." 

Another  approach  is  to  expand  components  of  the  vector  multipole  functions  in  the  a 
coordinate  system  into  expansions  containing  translation  coefficients  and  vector  spherical 
hannonics  in  the  p  coordinate  system.  The  expression  is  then  integrated  in  a  manner  similar  to  a 
Fourier  series  to  detennine  the  translation  coefficients.  The  fields  originating  from  sphere  a  are 
described  by  vector  multipole  functions  containing  terms  of  the  form  h(']( kra  )Y'nm  ( 0a ,  <pa  )  . 
Expanding  these  individual  components  in  the  P  coordinate  system  yields  the  following: 
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(23) 
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The  expansion  coefficients  (R a/})  are  functions  only  of  the  relative  position  vector 

Rap,  and  can  be  determined  from  the  following  integral: 

In  +1 

nS"  =  I  d(Pp  |  <i(cos  df}  )h{^ (kra  ) 

o  -i  (24) 


The  values  of  the  scattered  a/incident  p  multipole  fields  are  only  relevant  at  the  surface  of  sphere 
[!,  therefore  II  S”(R  op)  is  integrated  over  the  surface  of  sphere  [1  with  Op  and  cpp  as  our 
variables  of  integration. 

The  integral  is  solved  by  first  performing  the  dot  product  between  the  two  vector  spherical 
hannonics.  The  spherical  harmonic  terms  in  the  a  coordinate  system  are  then  expanded  in  terms 
of  bipolar  spherical  harmonics  containing  products  of  spherical  hannonics  in  the  [5  (local)  and  a,p 
(global)  coordinate  systems.  The  integration  is  then  performed  over  the  surface  of  sphere  P, 
and  the  summation  indices  are  simplified  with  the  use  of  orthogonality  conditions  imposed  by 
the  Clebsch-Gordan  coefficients.  The  final  form  for  the  expansion  coefficients  is  the  following: 
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(25) 


Since  the  expansion  coefficients  in  Eq.  (25)  only  translate  individual  multipole  terms  in  the 
U,  V,  and  W  wave  functions,  further  derivation  is  required  to  transform  the  entire  field 
expressions  as  given  by  Eqs.  (20)  and  (21).  To  detennine  S""1  and  T""' ,  it  is  first  necessary  to 
separate  the  spherical  Bessel  function  from  the  expansion  coefficients: 


K:z(Kap)=jAka«)z,c:(Kap) 


(26) 


This  provides  vector  multipole  tenns  of  the  form  j  ■  (ka  „  )Yvi  ( 6 p ,  (pp )  for  constructing  fields 

incident  to  sphere  P  (note  that  the  new  incident  fields  are  functions  of  spherical  Bessel  functions 
rather  than  spherical  Hankel  functions  of  the  second  kind): 
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(27) 
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Two  methods  can  then  be  used  to  derive  expressions  for  S™  and  T™‘ .  In  the  first  method, 
the  vector  fields  are  translated  directly  using  the  expansion  coefficients  of  Eq.  (27)  and  equating 
the  U,  V,  and  W  wave  functions  and  associated  scattered  wave  field  coefficients  with  the  U',  V', 
and  W'  wave  functions  and  incident  wave  field  coefficients.  In  the  second  method,  the  potentials 
for  the  fields  are  translated  using  Eq.  (27).  This  method  is  possible  since  the  fields  can  be 
derived  from  the  vector  potential  *F ,  and  requires  equating  the  potentials  and  corresponding 
amplitude  coefficients  and  subsequently  reconstructing  the  vector  fields  using  Eqs.  (8)  and  (9). 
Both  methods  produce  a  set  of  equivalent  solutions  for  the  translation  coefficients  S'™  and  T""‘ , 
the  simplest  of  which  are  the  following: 
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20 

Testing  of  Eqs.  (28-30)  showed  they  are  numerically  equivalent  to  those  of  Cruzan.“ 
Recurrence  formulas  have  been  derived  for  both  the  scalar  and  vector  addition  theorems,  and 

2 1  23 

significantly  reduce  the  number  of  calculations  for  computing  the  translation  coefficients.'  " 
However,  the  direct  expressions  in  Eqs.  (22)  and  (28-30)  were  used  to  demonstrate  proof-of- 
concept  for  the  modeling  approach. 


D.  Multiple-scattering  computations 


The  computations  are  performed  by  first  calculating  the  scattered  wave  fields  for  each 
particle  in  the  system  due  to  an  initial  plane  wave.  The  scattered  wave  fields  are  then  translated 
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between  all  of  the  particles  and  summed  for  each  particle.  A  new  incident  wave  field  (initial  + 
scattered  wave  fields)  is  then  used  to  compute  revised  scattered  wave  fields  from  each  particle. 
This  process  is  continued  by  iteration  until  the  scattered  wave  fields  converge  ( i.e .,  no  change  in 
amplitude  coefficients  between  consecutive  iterations).  Each  iteration  represents  a  successive 
order  of  scattering  (first  iteration  =  first-order  multiple  scattering,  second  iteration  =  second- 
order  multiple  scattering,  etc.).  Fig.  3  illustrates  the  computation  process. 


FIG.  3.  Flow  diagram  of  computation  steps  perfonned  in  the  elastic  wave  scattering  model. 

The  model  can  be  configured  to  provide  either  wave  field  images  or  spectral  results. 
Computation  limits  include  the  maximum  multipole  expansion  order  to  compute  (nmax),  and  the 
precision  limit  for  stopping  the  iterations.  Computer  algorithms  for  the  model  were  written, 
debugged,  and  compiled  in  Fortran  90.  Simulations  were  performed  on  a  personal  desktop 
computer  with  256  MB  RAM  and  a  1.7  GHz  processor  to  test  the  functioning  and  performance 
of  the  model.  Computation  times  varied  from  a  few  minutes  to  36  hours. 

In  addition  to  wave  field  images  and  attenuation  spectra,  it  is  also  possible  to  model  elastic 
wave  velocity  with  the  iterative  multipole  method.  Although  the  vector  multipole  functions  are 
time-independent  solutions  to  the  Navier  equation,  the  ability  to  manipulate  the  amplitude  of  the 
initial  plane  wave  as  a  function  of  frequency  allows  the  simulation  of  pulses  through  the  use  of 
Fourier  transforms.  This  is  accomplished  by  modifying  the  spectral  input  for  the  initial  plane 
wave  with  appropriate  wavelet  functions  [Fig.  4(a)],  calculating  the  interactions  of  the  modified 
wave  in  the  simulated  particle  pack,  and  applying  an  inverse  Fourier  transform  to  the  resultant 
spectra.  The  inverse  transform  creates  a  wave  pulse  in  the  time  domain.  Analysis  of  the  time 
delay  between  this  pulse  and  one  with  the  simulated  particle  pack  absent,  Fig.  4(b),  provides  the 
effective  wave  velocity  (note  that  pulse  attenuation  can  also  be  obtained).  This  procedure  is 
equivalent  to  simulating  the  frequency  characteristics  of  a  transmitting  ultrasonic  transducer. 

The  spatial  characteristics  of  the  transducer’s  wave  field  can  also  be  modeled  by  appropriately 
modifying  the  multipole  components  of  the  initial  wave.  For  example,  a  point  source  would 
produce  wave  fields  comprised  exclusively  of  monopole  and  dipole  moments  for  the  longitudinal 
and  shear  fields  respectively. 
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Time  (microseconds) 


(a) 


(b) 


FIG.  4.  Modulation  of  the  initial  plane  wave  spectrum  using  wavelet  transforms  (a),  and  time- 
delay  comparison  of  resultant  pulses  for  an  1 167-particle  dispersion  (b). 

III.  RESULTS 


A.  Simulations  of  two  to  sixteen  particles 

The  scattering  interactions  between  two  particles  were  initially  simulated  to  test  the  ability  of 
the  model  to  reproduce  wave  propagation  phenomena.  The  particles  were  modeled  as  quartz 
spheres  in  a  matrix  of  water  ice,  and  the  wave  field  patterns  were  imaged  to  reveal  whether 
realistic  wave  behavior  was  simulated.  The  acoustic  properties  of  quartz  and  ice  are  sufficiently 
different  to  clearly  show  scattering  phenomena  such  as  reflection,  but  not  too  different  to 
preclude  other  types  of  wave  propagation  such  as  refraction  through  particles.  Fig.  5  shows 
images  of  the  scattered  wave  fields  resulting  from  a  longitudinal  wave  incident  on  two  spheres 
aligned  along  the  direction  of  wave  propagation  (from  left  to  right  in  the  images).  The  incident 
wave  was  not  superimposed  on  the  scattered  wave  fields  in  order  to  emphasize  the  scattering 
mechanisms,  and  the  image  grey  scale  is  proportional  to  the  scattered  wave  field  amplitude.  The 
images  display  both  refractive  focusing  (forward  scattering)  of  the  longitudinal  wave,  Fig.  5(a), 
and  mode  conversion  of  a  small  part  of  the  longitudinal  wave  to  shear  waves,  Fig.  5(b). 
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(a)  (b) 

FIG.  5.  Wave  field  images  of  scattered  waves  resulting  from  an  incident  5-MHz  longitudinal 
wave  with  a  pair  of  1.0-mm  diameter  quartz  spheres  in  ice:  (a)  Scattered  longitudinal  wave 
displaying  focusing,  and  (b)  scattered  shear  wave  displaying  mode  conversion  from  longitudinal 
wave. 

Since  most  of  the  wave  behavior  in  Fig.  5  arises  from  single-sphere  scattering,  single¬ 
scattering  computations  were  compared  to  the  multiple-scattering  computations  to  ascertain  the 
differences  between  the  two  solutions  and  detennine  the  effects  of  multiple  scattering.  The 
single-scattering  computations  superimpose  the  single-scattering  solutions  from  each  of  the 
spheres  to  calculate  the  effective  field.  Figs.  6(a)  and  6(b)  are  difference  plots  of  the  multiply 
scattered  longitudinal  and  shear  waves  in  Figs.  5(a)  and  5(b),  respectively,  with  the  single¬ 
scattering  solutions.  The  difference  plots  highlight  the  spatial  regions  where  multiple  scattering 
has  the  largest  contribution  to  the  wave  fields.  The  images  in  Fig.  6  clearly  show  the  primary 
effect  of  multiple  scattering  in  the  two-sphere  configuration  is  the  shielding  of  the  second  (right) 
sphere  by  the  first  (left)  sphere. 
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(a)  (b) 

FIG.  6.  Displacement  differences  between  multiple-scattering  and  single-scattering  models  for 
an  incident  longitudinal  wave,  showing  a  scattered  longitudinal  wave,  z  component  (a),  and  a 
mode-converted  shear  wave,  x  component  (b). 

Simulations  were  subsequently  performed  for  2D  dispersions  containing  up  to  16  particles  to 
gauge  the  model’s  ability  to  calculate  multiple  scattering  in  more  complex  structures.  Fig.  7 
displays  the  propagation  of  a  plane  wave  with  mixed  longitudinal  and  shear  components  through 
an  ordered  [Fig.  7(a)]  and  random  [Fig.  7(b)]  configuration  of  quartz  particles.  Here  the  images 
show  the  x  component  of  the  shear  fields,  and  the  incident  wave  fields  are  superimposed  on  the 
scattered  wave  fields  to  give  a  truer  representation  of  the  elastic  fields  in  the  material  system. 
Figure  7(a)  shows  the  shear  wave  amplitude  decreasing  as  it  progresses  from  left  to  right  through 
the  particle  lattice.  A  total  decrease  in  peak-to-peak  amplitude  of  44%  is  observed  across  the 
image  field.  Since  the  constituent  material  properties  are  fully  elastic  with  no  built-in  attenuation 
properties,  the  attenuation  of  the  shear  waves  is  due  wholly  to  scattering.  The  attenuation 
probably  arises  from  an  acoustic  band  gap  since  the  waves  are  passing  though  an  ordered  array 
with  a  lattice  spacing  commensurate  with  the  wavelength.  In  contrast,  Figure  7(b)  shows 
localized  wave  field  concentrations  between  and  within  certain  particles.  This  localization  of 
wave  field  energy  may  be  due  to  resonance  effects.  Also,  the  wave  front  in  Fig.  7(b)  is  distorted 
from  the  initial  plane-wave  geometry  due  to  the  disordered  structure  of  the  particles,  in  contrast 
to  Fig.  7(a)  where  the  planar  wave  front  is  preserved. 
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(a)  (b) 

FIG.  7.  Total  (scattered  +  incident)  shear  wave  x  component  for  an  incident  mixed  longitudinal- 
shear  wave  at  1.0  MHz,  for  ordered  (a)  and  disordered  (b)  configurations  of  quartz  particles  in 
ice  (diameter  =  0.8- 1.4  mm). 

The  results  from  simulation  of  two  to  sixteen  particles  demonstrate  the  model’s  ability  to 
reproduce  basic  wave  propagation  behavior  such  as  refraction,  mode  conversion,  and  acoustic 
band  gaps  at  a  qualitative  level.  All  of  the  above  simulations  were  calculated  to  a  maximum 
multipole  order  of  nmax  =  6. 

B.  Simulations  of  dispersions  with  102-103  particles 

Spectra  and  wave  field  images  for  larger  3D  particle  packs  (100-700  particles)  were  also 
generated  to  demonstrate  the  ability  to  model  actual  material  microstructures.  Fig.  8  displays  the 
simulated  geometry  [Fig.  8(a)]  and  spectra  [Fig.  8(b)]  for  a  cylindrical  arrangement  of  200-pm 
diameter  NaCl  particles  in  a  soft  polymer  matrix.  For  the  simulation  of  spectra,  the  wave  fields 
are  evaluated  on  a  planar  array  of  points  with  a  circular  shape.  This  array  simulates  the  face  of 
an  ultrasonic  transducer  by  averaging  the  wave  field  amplitudes  over  a  defined  area.  This 
averaging  is  necessary  since  local  interference  effects  arise  when  the  fields  are  evaluated  at  a 
single  point.  The  particle  pack  dimensions  are  also  kept  constant  to  eliminate  effects  due  to 
changes  in  sample  size.  As  shown  in  the  spectra  [Fig.  8(b)],  the  model  predicts  that  the 
attenuation  increases  with  particle  volume  fraction  and  frequency.  These  results  are  both 
intuitive  (closer-packed  particles  lead  to  greater  ultrasonic  losses  due  to  scattering)  and  supported 
by  experiment. 
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Particle  pack: 

2.5-mm  diameter  x  1.25-mm  disk 
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FIG.  8.  Simulation  geometry  (a)  and  spectra  (b)  for  a  random  microstructure  of  200-pm 
diameter  NaCl  particles  in  a  soft  polymer,  displaying  attenuation  due  to  increased  volume 
fraction  of  particles. 


Fig.  9  compares  the  wave  fields  for  periodic  [Fig.  9(a)]  and  random  [Fig.  9(b)]  arrangements 
of  200-jum  diameter  particles  in  a  soft  polymer  matrix  at  50%  volume  fraction.  The  particle 
packings  are  three  dimensional  and  cylindrical  as  shown  in  Fig.  8(a),  with  the  image  plane 
slicing  through  the  center  of  the  cylinder  along  the  cylinder’s  axis.  The  incident  wave  is 
superimposed  in  these  images  in  order  to  evaluate  the  effect  of  order/disorder  on  the  total  elastic 
wave  field.  The  wave  field  images  in  Fig.  9  indicate  greater  fluctuation  and  non-unifonnity  for 
the  elastic  fields  associated  with  the  random  particle  configuration  [Fig.  9(b)]  as  compared  to  the 
periodic  lattice  [Fig.  9(a)],  The  ability  to  simulate  both  spectra  and  wave  field  patterns  allows  a 
more  thorough  characterization  of  wave  propagation  in  a  particle-filled  medium.  Both  spectral 
[Fig.  8(b)]  and  wave  field  simulations  [Fig.  9]  were  computed  to  a  maximum  multipole  order  of 

thnax  3 . 
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(a) 


(b) 


FIG.  9.  Wave  field  images  of  displacement  resulting  from  an  incident  0.5-MHz  longitudinal 
wave  in  a  cylindrical  packing  of  200-pm  diameter  NaCl  particles  in  rubber  at  50%  volume 
fraction:  (a)  Body-centered-cubic  (bee)  lattice  of  649  particles,  and  (b)  random  packing  of  736 
particles. 

The  ability  to  predict  elastic  wave  velocity  was  also  demonstrated  with  simulations  of 
dispersions  containing  600-pm  glass  particles  in  a  soft  polymer  matrix.  The  dispersions 
contained  up  to  3000  particles  and  50%  particle  volume  fraction.  The  results  show  an  initial 
increase  in  effective  velocity  with  increasing  packing  density  up  to  20%,  which  is  consistent  with 
experimental  observations.  However,  above  20%  particle  volume  the  model  under-predicts  the 
velocity.  These  results  demonstrate  the  viability  of  the  Fourier  transform  method  for  modeling 
wave  velocity  in  dispersions  of  spherical  particles,  but  also  indicate  the  current  simulations  are 
not  accurate  for  dense  packings.  This  deficiency  can  be  attributed  to  a  lack  of  convergence  for 
the  translational  addition  theorems  at  the  nmax  values  used  for  the  velocity  simulations  (nmax  =  4). 
The  next  section  examines  the  convergence  properties  of  the  model. 

C.  Convergence 

The  convergence  properties  of  the  model  can  be  divided  into  three  components: 

Convergence  of  the  addition  theorems,  convergence  of  the  spectra  of  the  effective  field,  and 
convergence  of  the  iterations.  These  components  behave  independent  of  one  another;  i.e.,  the 
convergence  of  one  component  is  not  dependent  on  the  convergence  of  the  other  two.  All  three 
of  these  convergences,  however,  are  dependent  on  the  maximum  computed  multipole  order 
(nmax),  on  the  frequency  of  the  propagating  wave,  and  on  the  structure  of  the  particle  dispersion 
(particle  distances,  particle  orientations,  and  degree  of  order-disorder).  Convergence  of  all  three 
components  is  necessary  for  computational  fidelity. 
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Addition  Theorem  Convergence 


The  convergence  of  the  addition  theorems  [Eqs.  (19-21)]  was  numerically  tested  by 
translating  the  elastic  fields  from  one  sphere  to  an  evaluation  point  on  another  sphere  and 
comparing  the  translated  fields  to  the  untranslated  fields  at  that  point.  Convergence  was 
evaluated  as  a  function  of  nmax,  frequency,  sphere  size,  sphere  separation,  orientation  of  the  two 
spheres,  and  orientation  of  the  evaluation  point  on  the  second  sphere.  Since  the  fields 
represented  elastic  fields  in  solids,  the  size  and  distance  scale  for  the  receiving  sphere  (sphere  P) 
was  in  the  centimeter  range,  and  the  frequencies  were  in  the  0-1.0  MHz  range. 

All  three  multipole  fields  (U,  V,  and  W)  were  tested  with  an  initial  quadrupole  (n= 2,  m= 1) 
moment.  A  quadrupole  moment  was  chosen  since  most  multiple-scattering  methods  seek  to 
model  many-body  wave  interactions  beyond  the  dipole  approximation,  which  becomes 
increasingly  less  accurate  as  scatterers  become  more  closely  spaced.  Additionally,  a  quadrupole 
moment  represents  a  good  compromise  between  n  >  2  moments  and  the  simple  dipole  moments 
since  computation  time  increases  by  approximately  the  sixth  power  of  the  multipole  order.  Also 
note  the  coupling  between  the  V  and  W  fields  in  the  addition  theorems  [Eqs.  (20-21)]  requires 
the  sum  of  the  two  transverse  fields  to  be  compared  instead  of  the  two  fields  individually  (i.e.,  V 
+  W  is  compared  to  V’  +  W’,  rather  than  V  to  V’  and  W  to  W’). 

Following  the  coordinate  conventions  in  Fig.  2,  Figs.  10-13  are  convergence  plots  for  sphere 
P  positioned  at  ©^p  =  34°,  ®ap  =  295°,  and  with  four  separation  distances  (f?ap  =  1-0,  2.0,  3.0, 
and  4.0  cm).  The  evaluation  point  on  sphere  P  was  positioned  at  0p  =  163°,  (pp  =  320°,  and  ap  = 
0.5  cm.  Each  figure  displays  results  for  a  different  frequency  (0.01,  0.10,  0.50,  and  1.00  MHz 
respectively).  At  low  frequencies  (Fig.  10),  the  results  show  that  field  translations  converge 
more  rapidly  as  a  function  of  nmax  for  larger  sphere  separations.  However,  as  the  field  frequency 
increases  to  0. 10  MHz,  the  convergence  curves  gradually  shift  and  begin  to  overlap  for  the  two 
most  distant  positions  (Fig.  1 1).  The  curves  continue  to  shift  with  higher  frequency  and  overlap 
for  all  positions  at  0.50  MHz  (Fig.  12).  The  curves  continue  to  shift,  however,  with  increasing 
frequency  and  cross  one  another  at  1.0  MHz  (Fig.  13). 
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0  5  10  15  20  25 

nmax 


FIG.  10.  Addition  theorem  convergence  for  translation  of  0.01 -MHz  shear  fields  between  two 
1 .0-cm  diameter  quartz  spheres  in  ice.  R  =  Rap  =  center-to-center  distance  between  spheres  in 
cm.  Numbers  in  parentheses  are  kRap  values. 


FIG.  1 1 .  Addition  theorem  convergence  for  translation  of  0. 10-MHz  shear  fields  between  two 
1 .0-cm  diameter  quartz  spheres  in  ice.  R  =  Rap  =  center-to-center  distance  between  spheres. 
Numbers  in  parentheses  are  kRap  values. 
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FIG.  12.  Addition  theorem  convergence  for  translation  of  0.50-MHz  shear  fields  between  two 
1 .0-cm  diameter  quartz  spheres  in  ice.  R  =  Rap  =  center-to-center  distance  between  spheres. 
Numbers  in  parentheses  are  kRap  values. 


FIG.  13.  Addition  theorem  convergence  for  translation  of  1.00-MHz  shear  fields  between  two 
1 .0-cm  diameter  quartz  spheres  in  ice.  R  =  Rap  =  center-to-center  distance  between  spheres. 
Numbers  in  parentheses  are  kRap  values. 
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The  convergence  curves  indicate  the  degree  of  error  present  in  the  field  translations  for  a 
specified  nmax,  frequency,  and  particle  separation.  For  example,  the  convergence  curves  for  Fig. 
12  would  correspond  to  the  two-sphere  simulations  in  Figs.  5  and  6  (since  the  frequency-sphere 
diameter  parameter  kd  is  equivalent).  From  Fig.  12,  however,  it  is  clear  computations  to  nmax  =  6 
will  result  in  field  translation  errors  on  the  order  of  10%  for  Rap  =  2,  and  an  nmax  value  of  12  is 
required  to  achieve  errors  on  the  order  of  1%.  Therefore,  although  the  iteration  sequence  of  the 
model  may  converge  to  a  solution,  the  field  translations,  which  are  crucial  to  accurately 
representing  the  multiple-scattering  interactions,  may  not  be  convergent  and  introduce  various 
degrees  of  error  into  the  solution. 

Spectral  Convergence 

In  multipole-based  optical  scattering  models  for  single  scatterers,  the  convergence  of  the 
scattered  wave  amplitudes  as  a  function  of  frequency  is  dependent  on  the  particle  diameter  d  and 
maximum  computed  multipole  order  nmax.  Simulation  results  verify  the  relationship  between 
nmax,  d,  and  convergence  frequency  vmax  (the  highest  frequency  at  which  the  wave  fields 
converge)  is  proportional  for  multiple-scattering  computations  of  elastic  waves  as  well: 


v 


max 


3  c 

- n 

An  d 


(31) 


The  symbol  c  is  the  wave  velocity  in  the  matrix.  An  example  is  the  following.  For  longitudinal 
wave  scattering  from  1-mm  diameter  quartz  particles  in  ice  (cv,  =  3. 98x10'’  cm/s)  and  for  nmax  = 

4,  the  maximum  frequency  at  which  the  spectrum  converges  is  vmax  ~  4  MHz.  This  result 
indicates  that  the  elastic  wave  spectra  will  be  convergent  for  frequencies  from  0-4  MHz  at  nmax  = 
4,  but  that  higher  nmax  values  will  be  required  for  convergence  at  higher  frequencies.  The 
spectral  convergence  is  independent  of  addition  theorem  and  iterative  convergence,  yet  lack  of 
spectral  convergence  can  contribute  significant  error  to  multiple-scattering  simulation  results. 

Iterative  Conversence 


Iterative  convergence  assures  that  the  back-and-forth  scattering  interactions  modeled  with 
each  iteration  step  converge  to  stable  wave  field  amplitudes.  This  is  physically  required  since 
each  successive  multiple  scattering  produces  a  diminishing  contribution  to  the  total  field, 
eventually  resulting  in  the  final  field  configuration.  Computationally,  however,  this  convergence 
is  not  always  assured  due  to  numerical  instabilities.  Simulations  of  a  variety  of  small  and  large 
particle  packings  have  shown  that  all  three  elastic  fields  converge  to  within  computational  limits 
(10'16  for  double  precision)  for  random  structures  and  a  wide  range  of  frequencies  [Fig.  14(a)]. 
Both  the  longitudinal  and  shear-electric  fields  also  converge  for  ordered  structures  such  as  cubic 
lattices  [Fig.  14(b)].  However,  the  shear  magnetic  field  does  not  converge  at  most  frequencies 
for  a  cubic  lattice,  but  rather  displays  fractional  differences  in  amplitude  on  the  order  of  10‘2  to 
10‘4  between  iteration  steps  [Fig.  14(b)].  The  origin  of  the  nonconvergence  for  the  shear- 
magnetic  field  has  not  yet  been  resolved.  However,  the  nonconvergence  has  not  affected  most  of 
the  multiple-scattering  simulations  to  date  due  to  its  small  effect  on  the  final  field  amplitudes  and 
its  limitation  to  ordered  lattices.  The  failure  of  the  shear-magnetic  field  to  converge  for  ordered 
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lattices  suggests  the  nonconvergence  may  be  attributable  to  resonance  effects,  numerical 
instabilities,  or  both. 


Iteration  Iteration 

(a)  (b) 

FIG.  14.  Iterative  convergence  of  simulations  for  a  cluster  of  eight  1.0-mm  diameter  quartz 
particles  in  ice  with  incident  2.5-MHz  longitudinal  and  shear  waves:  (a)  Random  cluster,  and  (b) 
cubic  cluster. 

In  addition  to  dispersion  structure  (particle-particle  orientations),  iterative  convergence  has 
also  been  found  to  be  dependent  on  frequency,  average  particle  separation  (volume  fraction  of 
particles),  and  nmax.  Iterative  convergence  appears  independent  of  both  addition  theorem  and 
spectral  convergence,  however,  indicating  that  all  three  convergence  criteria  must  be  considered 
in  iterative  multipole  models. 

IV.  DISCUSSION 

The  results  to  date  indicate  the  iterative  multipole  model  reproduces  realistic  wave 
propagation  behavior  at  a  qualitative  level.  Phenomena  such  as  mode  conversion,  refraction, 
scattering-induced  attenuation,  and  band-gap  attenuation  have  been  observed  in  both  wave  field 
images  and  spectra.  The  model  can  also  simulate  large  numbers  of  particles  and  complex 
random  structures  directly.  Random  particle  packings  of  up  to  3000  particles  and  50%  particle 
volume  fraction  have  been  modeled.  The  iterations  are  generally  numerically  stable  and 
converge  for  most  particle  configurations  and  fields.  The  one  exception  that  has  been  found  is 
for  the  shear-magnetic  field  in  cubic  lattices,  where  the  convergence  of  the  iterations  plateaus  at 
an  error  of  approximately  10" . 

Addition  theorem  convergence  is  the  current  limiting  factor  for  the  iterative  multipole 
approach  implemented  on  a  personal  computer.  Due  to  the  intensive  computations  involved  with 
the  addition  theorems,  the  maximum  computed  multipole  order  (nmax)  has  so  far  been  limited  to 
values  of  3-6.  Convergence  studies  of  the  addition  theorems  and  model  comparisons  with 
experimental  data  indicate  that  these  values  are  sufficient  for  dilute  particle  dispersions  (up  to 
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20%  particle  volume  fraction),  but  are  too  low  for  denser  packings.  Consequently,  simulations 
performed  to  date  with  high  particle  densities  probably  have  significant  error  in  the  multiple¬ 
scattering  contributions  due  to  inadequate  convergence  of  the  field  translations.  The  numerical 
stability  of  the  iterative  procedure  is  surprisingly  robust,  however,  with  iterative  convergence 
unaffected  by  either  addition  theorem  convergence  or  spectral  convergence.  Care  must  be  taken, 
therefore,  in  interpreting  the  accuracy  of  the  model  results  with  regards  to  iterative  convergence. 

A  common  approach  to  validating  numerical-based  multiple-scattering  models  is  to  compare 
the  results  for  a  small  cluster  of  spheres  to  other  numerical  approaches  such  as  the  boundary 
element  method.8  This  approach  is  feasible  for  clusters  of  two  or  three  spheres,  but  often 
becomes  computationally  intractable  for  large  particle  packings  with  more  than  a  few  dozen 
particles.  Alternatively,  model  results  can  be  compared  to  experimental  data.  This  approach  has 
been  used  extensively  for  evaluating  classical  formulations  where  multiple  scattering  is  treated 
approximately.32'33  Analogous  experimental  validation  of  the  iterative  multipole  method  is 
therefore  an  appropriate  test  for  the  model. 

The  current  limitations  of  the  iterative  multipole  model  restrict  its  validity  to  dilute 
dispersions.  Here  the  particles  are  spaced  sufficiently  to  allow  the  use  of  low  nmax  values  for 
addition  theorem  convergence.  However,  it  is  difficult  to  assess  the  accuracy  of  multiple¬ 
scattering  computations  for  these  conditions  since  single  scattering  dominates  in  dilute 
dispersions.  One  useful  tool  in  evaluating  the  accuracy  of  the  multiple-scattering  computations 
is  to  compare  their  results  to  those  of  single-scattering  computations  for  the  same  particle 
configuration.  Such  comparisons  show  that  the  iterative  multipole  model  agrees  with  the  single¬ 
scattering  model  for  low  nmax  values.  This  is  a  consistent  finding  for  low  particle  densities  since 
the  multiple-scattering  contributions  to  the  total  elastic  wave  field  will  be  small  in  comparison  to 
the  single-scattering  contributions.  It  is  also  a  consistent  finding  for  high  particle  densities  since 
it  shows  the  low  nmax  values  are  insufficient  for  computing  accurate  translation  coefficients  for 
multiple  scattering,  again  rendering  their  contribution  to  the  total  elastic  wave  field  negligible. 

The  difficulty  in  validating  the  iterative  multipole  model  with  dilute  particle  packings 
extends  to  other  multiple-scattering  models  as  well,  and  additionally  illustrates  how  single 
scattering  can  overwhelm  contributions  from  multiple  scattering.  In  such  cases  multiple¬ 
scattering  models  may  produce  accurate  results  for  a  particle  dispersion,  although  the  multiple 
scattering  may  be  insufficiently  or  incorrectly  calculated.  Validation  with  dense  particle  systems 
is  therefore  necessary  for  accurately  evaluating  a  model’s  performance. 

The  validation  results  to  date  for  the  iterative  multipole  model  indicate  that  although  the 
model  appears  to  be  running  properly  and  accurately  for  dispersions  with  low  particle  densities, 
the  computational  limitation  on  nmax  renders  the  model  quantitatively  inaccurate  for  dispersions 
with  high  particle  densities.  Fig.  15  illustrates  this  conclusion  by  comparing  longitudinal  wave 
velocity  results  for  the  iterative  multipole  model,  the  single-scattering  model,  and  experimental 
measurements  for  600-pm  glass  beads  in  an  elastic  polymer  matrix.  Note  the  agreement  between 
the  iterative  multipole  model  and  single-scattering  model  at  high  particle  densities,  indicating  the 
insufficient  computation  of  multiple  scattering  due  to  low  nmax. 
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FIG.  15.  Comparison  of  longitudinal  wave  velocity  at  0.5  MHz  for  the  iterative  multipole 
model,  the  single-scattering  model,  and  experimental  measurements  for  600-pm  glass  beads  in 
an  elastic  polymer  matrix. 

Several  solutions  are  available  to  resolve  the  computational  limits  to  the  iterative  multipole 
method.  These  include  increasing  the  speed  and  efficiency  of  the  addition  theorem  computations 
with  the  use  of  recurrence  algorithms.  '  The  computations  can  also  be  simplified  for  fluid  and 
mechanically  soft  matrices  by  limiting  the  propagating  elastic  fields  to  longitudinal  waves. 
Finally,  particle-particle  interactions  that  do  not  significantly  contribute  to  the  multiple-scattering 
process  (for  example,  from  widely  separated  particles)  can  be  omitted,  thereby  reducing  the 
number  of  field  translation  computations.  Other  foreseeable  improvements  for  the  iterative 
multipole  approach  consist  of  expanding  the  ability  to  model  more  complex  particulate  systems. 
These  include  simulating  viscoelastic  media  with  the  use  of  complex  material  properties  and 
wave  vectors,  and  simulating  mixtures  of  spherical  voids  and  particles. 

In  comparison  to  other  multipole  approaches,  the  iterative  or  order-of-scattering  method 
appears  to  be  a  more  direct  and  accessible  approach  for  modeling  elastic  waves  scattering  in 
arbitrary  dispersions  of  spherical  particles.  The  iterative  multipole  method  is  numerically  stable, 
particularly  for  random  particle  configurations,  whereas  stability  cannot  be  guaranteed  for 
inversion  of  large  matrices.  Finally,  the  iterative  multipole  method  is  a  highly  generalized 
approach  that  can  be  applied  to  dispersions  of  spherical  scatterers  of  almost  any  configuration 
and  composition.  No  fundamental  assumptions  or  approximations  about  the  dispersion  structure 
are  made.  For  example,  the  dispersion  structure  is  not  assumed  to  be  random  or  homogeneous  at 
large  scale  sizes  as  required  by  fast  multipole  methods.  Additionally,  the  structure  does  not  have 
to  be  built  of  repeating  units  or  limited  to  two  dimensions  as  is  often  the  practice  for  finite 
methods. 
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V.  CONCLUSIONS 


An  iterative  multipole  method  has  been  developed  to  model  arbitrary  configurations  of 
spherical,  elastic  particles  in  an  elastic  matrix,  and  has  been  implemented  on  a  personal  computer 
to  test  its  functioning  and  performance.  Initial  results  for  systems  of  two  to  several  hundred 
particles  show  that  the  model  predicts  many  aspects  of  elastic  wave  behavior  at  the  qualitative 
level.  The  model  is  also  quantitatively  accurate  for  dispersions  with  low  particle  densities 
(<20%  particle  volume  fraction).  However,  assessing  the  accuracy  of  a  multiple-scattering 
model  for  dilute  dispersions  is  difficult  since  single-scattering  interactions  dominate.  At  higher 
particle  densities,  the  model  predictions  deviate  due  to  limitations  in  computing  sufficiently  high 
multipole  orders  (nmax  =  8-16)  for  convergence  of  the  addition  theorems.  Methods  for  resolving 
these  limitations  include  the  use  of  recurrence  formulas  for  computing  the  addition  theorems, 
approximations  that  reduce  the  number  of  computations  by  limiting  the  multiple  scattering  to  the 
most  significant  interactions,  and  implementation  of  the  model  on  larger  computer  systems. 
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APPENDIX 

The  //  functions  for  solution  of  the  elastic  wave  boundary  condition  matrix,  Eq.  (18)  are 
presented  below.  The  // ’s  are  functions  of  the  multipole  order  N,  longitudinal  wave  vector  kL, 
shear  wave  vector  ks,  and  the  sphere  radius  a.  The  spherical  radial  functions  are  denoted  by  z, 
and  vary  according  to  the  radial  function  designated  in  Eq.  (18)  [j  =  jn  (kr) ,  h  =  ( kr ) ,  and 

g  =  h^(kr)]. 


=  N  \  L  ZN+l(kLa)> 

(Al) 

kLci 

r/2(z) 

=  iJ(N)(N  +  \)  Z^ksCl)  , 
ksa 

(A2) 

*h(z) 

=  zN(ksa), 

(A3) 

V4(z) 

_  ZN(kLa ) 

k,  a 

(A4) 

Approved  for  public  release;  distribution  unlimited. 


Iterative  simulation  of  scattering  26 


I 


r/5(z)  = 


J(N)(N  + 1) 


(A  +  l) 


zN(ksa) 

ksa 


-  z 


N+ 1 


(ksa) 


(A5) 


76(z)  =  ~^kLz  N(kLa) 

(N-l)(N)  N  +  l 


+  2  /ukL 


(, kLa )2  N  +  2 


zN(kLa)  + 


zN+i(kLa )  , 


kLa 


T) 


(z)  =  2i/jks  J(N)(N  + 1) 


(A-l) 


zN(ksa)  z 

N+ 1  (V*) 

(ksa)2  ksa 


r/s(z)  =  juks 


zN(ksa) 


(N-l)^^~ZN+l(ksa) 

ksa 


tj9(z)  =  2juk, 


(N- 1) 


zN(kLa) 

_  ZN+\  (kLa) 

(kLa)2  kLa 


n  io(z)  = 


ijuks 


J(N)(N  + 1) 

2(A^2 -1)  TV  + 1 


{ksaY 


N  +  2 


'N 


(ksa)  + 


'N+ 1 


(M) 


ksa 


(A6) 

(A7) 

(A8) 

(A9) 

(A10) 


REFERENCES 

1.  R.  Marklein,  K.  J.  Langenberg,  R.  Barmann,  K.  Mayer,  and  S.  Klaholz,  “Nondestructive 
Testing  with  Ultrasound:  Numerical  Modeling  and  Imaging,”  in  Acoustical  Imaging,  Vol. 
22,  edited  by  P.  Tortoli  and  L.  Masotti  (Plenum  Press,  New  York,  1996),  pp.  757-764. 

2.  A.  A.  Gusev,  “Numerical  identification  of  the  potential  of  whisker  and  platelet  filled 
polymers,”  Macromolecules  34,  3081-3093  (2001). 

3.  V.  E.  Ostashev,  D.  K.  Wilson,  L.  Liu,  D.  F.  Aldridge,  N.  P.  Symons,  and  D.  Marlin, 
“Equations  for  finite-difference,  time  domain  simulation  of  sound  propagation  in  moving 
inhomogeneous  media  and  numerical  implementation,”  J.  Acoust.  Soc.  Am.  117,  503-517 
(2005). 


Approved  for  public  release;  distribution  unlimited. 


Iterative  simulation  of  scattering  27 


4.  A-K.  Hamid,  I.  R.  Ciric,  and  M.  Hamid,  “Iterative  solution  of  the  scattering  by  an  arbitrary 
configuration  of  conducting  or  dielectric  spheres,”  IEE  Proceedings-H  138,  565-572  (1991). 

5.  K.  A.  Fuller  and  G.  W.  Kattawar,  “Consummate  solution  to  the  problem  of  classical 
electromagnetic  scattering  by  an  ensemble  of  spheres.  II.  Clusters  of  arbitrary  configuration,” 
Opt.  Lett.  13,  1063-1065  (1988). 

6.  D.  W.  Mackowski,  “Analysis  of  radiative  scattering  for  multiple  sphere  configurations,” 

Proc.  R.  Soc.  Lond.  A  433,  599-614  (1991). 

7.  S.  Koc  and  W.  C.  Chew,  “Calculation  of  acoustical  scattering  from  a  cluster  of  scatterers,”  J. 
Acoust.  Soc.  Am.  103,  721-734  (1998). 

8.  N.  A.  Gumerov  and  R.  Duraiswami,  “Computation  of  scattering  from  N  spheres  using 
multipole  reexpansion,”  J.  Acoust.  Soc.  Am.  112,  2688-2701  (2002). 

9.  Z.  Liu,  C.  T.  Chan,  P.  Sheng,  A.  L.  Goertzen,  and  J.  H.  Page,  “Elastic  wave  scattering  by 
periodic  structures  of  spherical  objects:  Theory  and  experiment,”  Phys.  Rev.  B  62,  2446- 
2457  (2000). 

10.  L.  Greengard  and  V.  Rokhlin,  “A  fast  algorithm  for  particle  simulations,”  J.  Comput.  Phys. 
73,325-348  (1987). 

1 1.  W.  Greiner  and  J.  A.  Maruhn,  Nuclear  Models  (Springer- Verlag,  Berlin,  1996). 

12.  K.  S.  Thorne,  “Multipole  expansions  of  gravitational  radiation,”  Rev.  Mod.  Phys.  52,  299- 
339  (1980). 

13.  D.  A.  Varshalovich,  A.  N.  Moskalev,  and  V.  K.  Khersonskii,  Quantum  Theory  of  Angular 
Momentum,  (World  Scientific,  Singapore,  1988). 

14.  J.  A.  Stratton,  Electromagnetic  Theory  (McGraw-Hill,  New  York,  1941). 

15.  M.  E.  Rose,  Elementary  Theory  of  Angular  Momentum  (John  Wiley  and  Sons,  New  York, 
1957). 

16.  A.  R.  Edmonds,  Angular  Momentum  in  Quantum  Mechanics  (Princeton  University  Press, 
Princeton,  New  Jersey,  1957). 

17.  J.  D.  Jackson,  Classical  Electrodynamics,  second  edition  (John  Wiley  and  Sons,  New  York, 
1975). 

18.  C.  F.  Ying  and  R.  Truell,  “Scattering  of  a  plane  longitudinal  wave  by  a  spherical  obstacle  in 
an  isotropically  elastic  solid,”  J.  Appl.  Phys.  27,  1086-1097  (1956). 


Approved  for  public  release;  distribution  unlimited. 


Iterative  simulation  of  scattering  28 


19.  L.  Knopoff,  “Scattering  of  shear  waves  by  spherical  obstacles,”  Geophysics  24,  209-219 
(1959). 

20.  L.  Knopoff,  “Scattering  of  compression  waves  by  spherical  obstacles,”  Geophysics,  24,  SO¬ 
SO  (1959). 

21.  N.  G.  Einspruch,  E.  J.  Witterholt,  and  R.  Truell,  “Scattering  of  a  plane  transverse  wave  by  a 
spherical  obstacle  in  an  elastic  medium,”  J.  Appl.  Phys.  31,  806-818  (1960). 

22.  D.  L.  Jain  and  R.  P.  Kanwal,  “Scattering  of  P  and  S  waves  by  spherical  inclusions  and 
cavities,”  J.  Sound  Vib.  57,  171-202  (1978). 

23.  M.  K.  Hinders,  “Plane-elastic-wave  scattering  from  an  elastic  sphere,”  Nuovo  Cimento  B 
106B,  799-817(1991). 

24.  B.  Friedman  and  J.  Russek,  “Addition  theorems  for  spherical  waves,”  Q.  Appl.  Math.  12,  13- 
23  (1954). 

25.  M.  E.  Rose,  “The  electrostatic  interaction  of  two  arbitrary  charge  distributions,”  J.  Math. 
Phys.  (Cambridge,  Mass.)  37,  215-222  (1958). 

26.  S.  Stein,  “Addition  theorems  for  spherical  wave  functions,”  Q.  Appl.  Math.  19,  15-24  (1961). 

27.  R.  A.  Sack,  “Three-dimensional  addition  theorem  for  arbitrary  functions  involving 
expansions  in  spherical  harmonics,”  J.  Math.  Phys.  (N.Y.)  5,  252-259  (1964). 

28.  O.  R.  Cruzan,  “Translational  addition  theorems  for  spherical  vector  wave  functions,”  Q. 

Appl.  Math.  20,  33-40(1962). 

29.  W.  C.  Chew,  “Recurrence  relations  for  three-dimensional  scalar  vector  addition  theorem,”  J. 
Electromagn.  Waves  Appl.  6,  133-142  (1992). 

30.  W.  C.  Chew  and  Y.  M.  Wang,  “Efficient  ways  to  compute  the  vector  addition  theorem,”  J. 
Electromagn.  Waves  Appl.  7,  651-665  (1993). 

31.  D.  W.  Mackowski,  “Calculation  of  total  cross  sections  of  multiple-sphere  clusters,”  J.  Opt. 
Soc.  Am.  A  11,  2851-2861  (1994). 

32.  L.  W.  Anson  and  R.  C.  Chivers,  “Ultrasonic  velocity  in  suspensions  of  solids  in  solids — a 
comparison  of  theory  and  experiment,”  J.  Phys.  D:  Appl.  Phys.  26,  1506-1575  (1993). 

33.  R.  E.  Challis,  J.  S.  Tebbutt,  and  A.  K.  Holmes,  “Equivalence  between  three  scattering 
formulations  for  ultrasonic  wave  propagation  in  particulate  mixtures,”  J.  Phys.  D:  Appl. 
Phys.  31,  3481-3497  (1998) 


Approved  for  public  release;  distribution  unlimited. 


Iterative  simulation  of  scattering  29 


